#---------------------------------------------------------------------------
# Statistical results in appendix of "Bypass Aid and Unrest in Autocracies"
# Matthew DiLorenzo
# mdilorenzo@wm.edu
# Last updated 08-21-2017
#---------------------------------------------------------------------------

## Load packages
require(foreign)
require(MASS)
require(stargazer)
require(lmtest)
require(sandwich)
require(car)
require(dplyr)
require(erer)

## Load two custom functions: empiricalMarginsNB() and twostageboot().
source("~/Dropbox/article-manuscripts/baua-isq/baua-isq-revised/scripts/baua-custom-functions.R")

## Load in data set
load(file="/Users/matthewdilorenzo/Dropbox/Dissertation/Data/baas-data-10-21-2016.RData")


## Correcting for scale of variables
dat$bypass_pop_log <- log(
  ((dat$bypass2_t_t1*1000000) / (dat$population_t1*1000)) + 1
  )
dat$public_pop_log <- log(
  ((dat$public_t_t1*1000000) / (dat$population_t1*1000)) + 1
  )



#---------------------------------------------------------------
# Measure of unrest excluding riots
#---------------------------------------------------------------
dat$no_riots<-dat$banks_domestic2 + dat$banks_domestic8
dat$no_riots_dummy <- ifelse(dat$no_riots > 0, 1, 0)

## Model 1
m1 <- glm.nb(no_riots ~ bypass_pop_log +
               gov_index_t1 + 
               civil_conflict_t1 + 
               nd_count_t1,
             data=dat, 
             na.action="na.omit")

## Model 2
m2 <- update(m1, . ~ pct_bypass2_p_t1 + . - bypass_pop_log)

## Model 3
m3 <- update(m1, . ~ . + factor(cowcode) + factor(year) - 1)

## Model 4
m4 <- update(m3, . ~ . + public_pop_log)

## Model 5
m5 <- update(m2, . ~ . + factor(cowcode) + factor(year) - 1)
summary(m5)
coeftest(m5, vcovHC(m5, "HC1"))


## Model 6
lpm1 <- lm(no_riots_dummy ~ bypass_pop_log + 
             gov_index_t1 + 
             civil_conflict_t1 + 
             nd_count_t1 +
             public_pop_log + 
             factor(cowcode) + 
             factor(year) - 1,
           data=dat, 
           na.action="na.omit")
summary(lpm1)
coeftest(lpm1, vcovHC(lpm1, "HC1"))
## Model 7
lpm2 <- update(lpm1, . ~ pct_bypass2_p_t1 + . - bypass_pop_log)

## Collect robust standard errors in list for table
models <- list(m1, m2, m3, m4, m5, lpm1, lpm2)
ses <- lapply(models, function(x) sqrt(diag(vcovHC(x, type = "HC1"))))

## Variable names for covariates
covariates <- c("Bypass Aid / Capita",
                "Bypass Ratio",
                "Governance Index",
                "Civil Conflicts",
                "Natural Disasters",
                "Government Aid / Capita")

## Export to table format 
sink("~/Dropbox/article-manuscripts/baua-isq/baua-isq-final/baua-table-appx-riots.tex")
stargazer(m1, m2, m3, m4, m5, lpm1, lpm2,
          se = ses,
          label = "riot-models",
          title = "Bypass Aid and Domestic Unrest (Excluding Riots), 2005-2010",
          column.labels = c("\\emph{Model 1}", "\\emph{Model 2}", 
                            "\\emph{Model 3}", "\\emph{Model 4}", 
                            "\\emph{Model 5}", "\\emph{Model 6}",
                            "\\emph{Model 7}"),
          dep.var.caption = c("Unrest Events"),
          dep.var.labels.include = F,
          covariate.labels = covariates,
          omit = c("cowcode", "year", "Constant"),
          p.auto = TRUE,
          t.auto = TRUE,
          digits = 2,
          column.sep.width = c("0pt"),
          add.lines = list(c("Country-fixed Effects", rep("N", 2) ,rep("Y", 5)),
                           c("Year-fixed Effects", rep("N", 2) ,rep("Y", 5))),
          notes.append = TRUE,
          notes.align = "l",
          notes.label = "",
          font.size = "footnotesize",
          notes = c("Two-tailed tests. (\\possessivecite{White1980} HC1 standard errors)", 
                    "Negative binomial regression (Models 1-5) and OLS (Models 6 and 7). Outcome in Models 6 and 7 is dichotomous indicator of at least one unrest event."),
          float.env = "sidewaystable",
          model.names = FALSE)
sink()






#---------------------------------------------------------------
# Logit models
#---------------------------------------------------------------

## Model 1
m1 <- glm(unrest_dummy ~ bypass_pop_log +
            gov_index_t1 + 
            civil_conflict_t1 + 
            nd_count_t1,
          family = binomial(link = "logit"),
          data=dat, 
          na.action="na.omit")

## Model 2
m2 <- update(m1, . ~ pct_bypass2_p_t1 + . - bypass_pop_log)

## Model 3
m3 <- update(m1, . ~ . + factor(cowcode) + factor(year) - 1)

## Model 4
m4 <- update(m3, . ~ . + public_pop_log)

## Model 5
m5 <- update(m2, . ~ . + factor(cowcode) + factor(year) - 1)
summary(m5)


## Variable names for covariates
covariates <- c("Bypass Aid / Capita",
                "Bypass Ratio",
                "Governance Index",
                "Civil Conflicts",
                "Natural Disasters",
                "Government Aid / Capita")

## Export to table format 
sink("~/Dropbox/article-manuscripts/baua-isq/baua-isq-final/baua-table-appx-logit.tex")
stargazer(m1, m2, m3, m4, m5,
          label = "logit-models",
          title = "Bypass Aid and Domestic Unrest (Logit Models), 2005-2010",
          column.labels = c("\\emph{Model 1}", "\\emph{Model 2}", 
                            "\\emph{Model 3}", "\\emph{Model 4}", 
                            "\\emph{Model 5}"),
          dep.var.caption = c("Unrest Events"),
          dep.var.labels.include = F,
          covariate.labels = covariates,
          omit = c("cowcode", "year", "Constant"),
          p.auto = TRUE,
          t.auto = TRUE,
          digits = 2,
          column.sep.width = c("0pt"),
          add.lines = list(c("Country-fixed Effects", rep("N", 2) ,rep("Y", 3)),
                           c("Year-fixed Effects", rep("N", 2) ,rep("Y", 3))),
          notes.append = TRUE,
          notes.align = "l",
          notes.label = "",
          font.size = "footnotesize",
          notes = c("Two-tailed tests. Estimated standard errors in parentheses.", 
                    "Logit regression (Models 1-5). Outcome is dichotomous indicator of at least one unrest event."),
          float.env = "sidewaystable",
          model.names = FALSE)
sink()






#---------------------------------------------------------------
# Zero-inflated negative binomial regression models
#---------------------------------------------------------------
require(ggplot2)
require(pscl)
require(MASS)
require(boot)


zdat <- dat %>%
  dplyr::select(unrest, 
                bypass_pop_log, 
                pct_bypass2_p_t1, 
                gov_index_t1,
                civil_conflict_t1,
                nd_count_t1,
                year,
                cowcode) %>%
  na.exclude()


m1.zinb <- zeroinfl(unrest ~ I(bypass_pop_log*100) +
                      gov_index_t1 + 
                      civil_conflict_t1 + 
                      nd_count_t1 + 
                      factor(year) | 
                      factor(cowcode), 
                    data = zdat, 
                    dist = "negbin",
                    na.action = "na.omit")
summary(m1.zinb)



m2.zinb <- zeroinfl(unrest ~ pct_bypass2_p_t1 +
                      gov_index_t1 + 
                      civil_conflict_t1 + 
                      nd_count_t1 +
                      factor(year) | 
                      factor(cowcode), 
                    data = zdat, 
                    dist = "negbin",
                    na.action = "na.omit")
summary(m2.zinb)


## Variable names for covariates
covariates <- c("Bypass Aid / Capita",
                "Bypass Ratio",
                "Governance Index",
                "Civil Conflicts",
                "Natural Disasters")

## Export to table format 
sink("~/Dropbox/article-manuscripts/baua-isq/baua-isq-final/baua-table-appx-zinb.tex")
stargazer(m1.zinb, m2.zinb,
          label = "zinb-models",
          title="Zero-inflated Negative Binomial Regression Models",
          column.labels = c("\\emph{Model 1}", "\\emph{Model 2}"),
          dep.var.caption = c("Unrest Events"),
          dep.var.labels.include = F,
          covariate.labels = covariates,
          omit = c("cowcode", "year", "Constant"),
          p.auto = TRUE,
          t.auto = TRUE,
          digits = 2,
          column.sep.width = c("0pt"),
          add.lines = list(c("Year-fixed Effects", rep("Y", 2))),
          notes.append = TRUE,
          notes.align = "l",
          notes.label = "",
          font.size = "footnotesize",
          notes = c("Two-tailed tests. Estimated standard errors in parentheses.", 
                    "Zero-inflated negative binomial regression."),
          model.names = FALSE)
sink()








#---------------------------------------------------------------
# Changes in military expenditures
#---------------------------------------------------------------
dat$milex_change <- dat$milex - dat$milex_lag

milex1 <- lm(milex_change ~ bypass_pop_log + 
               interstate_crises_t1 + 
               civil_conflict_t1 + 
               factor(year),
             data = dat)
milex1.hcses <- coeftest(milex1, vcov = vcovHC(milex1, type = "HC1"))
milex1.hcses



milex2 <- lm(milex_change ~ pct_bypass2_p_t1 + 
               interstate_crises_t1 + 
               civil_conflict_t1 + 
               factor(year),
             data = dat)
milex2.hcses <- coeftest(milex2, vcov = vcovHC(milex2, type = "HC1"))
milex2.hcses


covariates <- c("Bypass Aid / Capita",
                "Bypass Ratio",
                "International Disputes",
                "Civil Conflicts")

## Export to table format 
sink("~/Dropbox/article-manuscripts/baua-isq/baua-isq-final/baua-table-appx-milex.tex")
stargazer(milex1.hcses, milex2.hcses,
          label = "milex-models",
          title="Bypass Aid and Change in Military Expenditures (\\% Gov. Spending)",
          column.labels = c("\\emph{Model 1}", "\\emph{Model 2}"),
          dep.var.caption = c("Dependent variable: Change in Military Expenditures"),
          dep.var.labels.include = F,
          covariate.labels = covariates,
          omit = c("cowcode", "year", "Constant"),
          p.auto = TRUE,
          t.auto = TRUE,
          digits = 2,
          column.sep.width = c("0pt"),
          add.lines = list(c("Year-fixed Effects", rep("Y", 2))),
          notes.append = TRUE,
          notes.align = "l",
          notes.label = "",
          font.size = "footnotesize",
          notes = c("Two-tailed tests. Estimated standard errors in parentheses.", 
                    "OLS estimates."),
          model.names = FALSE)
sink()





#---------------------------------------------------------------
# AidData Models
#---------------------------------------------------------------

## Create bypass/government aid per capita variables 
dat$ad_bypass_pop_log <- log(dat$ad_bypass_total_t1/(dat$population_t1*1000) + 1)
dat$ad_public_pop <- dat$ad_govt_total_t1/dat$population_t1

## Model 1
a1 <- glm.nb(unrest ~ ad_bypass_pop_log +
               gov_index_t1 + 
               civil_conflict_t1 + 
               nd_count_t1,
             data=dat, 
             na.action="na.omit")

## Model 2
a2 <- update(a1, . ~ ad_bypass_proportion_t1 + . - ad_bypass_pop_log)


## Collect robust standard errors in list for table
models <- list(a1, a2)
ses <- lapply(models, function(x) sqrt(diag(vcovHC(x, type = "HC1"))))

covariates <- c("Bypass Aid / Capita",
                "Bypass Ratio",
                "Governance Index",
                "Civil Conflicts",
                "Natural Disasters")

## Export to table format 
sink("~/Dropbox/article-manuscripts/baua-isq/baua-isq-final/baua-table-appx-aiddata.tex")
stargazer(a1, a2,
          se = ses,
          label = "ad-models",
          title="Alternative AidData Measure",
          column.labels = c("\\emph{Model 1}", "\\emph{Model 2}"),
          dep.var.caption = c("Dependent variable: Unrest Events"),
          dep.var.labels.include = F,
          covariate.labels = covariates,
          omit = c("cowcode", "year", "Constant"),
          p.auto = TRUE,
          t.auto = TRUE,
          digits = 2,
          column.sep.width = c("0pt"),
          notes.append = TRUE,
          notes.align = "l",
          notes.label = "",
          font.size = "footnotesize",
          notes = c("Two-tailed tests. (\\possessivecite{White1980} HC1 standard errors)", 
                    "Negative binomial regression models."),
          model.names = FALSE)
sink()





#---------------------------------------------------------------
# Alternative operationalizations of bypass aid
#---------------------------------------------------------------

## Model 1
ma1<-glm.nb(unrest ~ pct_bypass1_p_t1 +
              gov_index_t1 +
              civil_conflict_t1 +
              nd_count_t1 +
              factor(cowcode) +
              factor(year), data=dat, na.action="na.omit")

## Model 2
ma2 <- update(ma1, . ~ bypass_pct_gdp + . - pct_bypass1_p_t1)


## Collect robust standard errors in list for table
models <- list(ma1, ma2)
ses <- lapply(models, function(x) sqrt(diag(vcovHC(x, type = "HC1"))))

## Variable names for covariates
covariates <- c("Bypass Ratio (Excluding Public-Private Partnerships)",
                "Bypass / GDP",
                "Governance Index",
                "Civil Conflicts",
                "Natural Disasters")

## Export to table format 
sink("~/Dropbox/article-manuscripts/baua-isq/baua-isq-final/baua-table-appx-alt-oecd.tex")
stargazer(ma1, ma2,
          se = ses,
          label = "alt-oecd-models",
          title="Alternative OECD Measures of Bypass Aid",
          column.labels = c("\\emph{Model 1}", "\\emph{Model 2}"),
          dep.var.caption = c("Unrest Events"),
          dep.var.labels.include = F,
          covariate.labels = covariates,
          omit = c("cowcode", "year", "Constant"),
          p.auto = TRUE,
          t.auto = TRUE,
          digits = 2,
          column.sep.width = c("0pt"),
          add.lines = list(c("Country-fixed Effects", rep("Y", 2)), c("Year-fixed Effects", rep("Y", 2))),
          notes.append = TRUE,
          notes.align = "l",
          notes.label = "",
          font.size = "footnotesize",
          notes = c("Two-tailed tests. (\\possessivecite{White1980} HC1 standard errors)", 
                    "Negative binomial regression."),
          model.names = FALSE)
sink()


#---------------------------------------------------------------
# Controlling for Strategic Importance of Recipient Country
#---------------------------------------------------------------
names(dat)
## Model 1
m1 <- glm.nb(unrest ~ bypass_pop_log +
               gov_index_t1 + 
               civil_conflict_t1 + 
               nd_count_t1 +
               yrs_since_independence +
               us_ally +
               n_maj_allies +
               opec_member,
             data=dat, 
             na.action="na.omit", maxit = 100)
summary(m1)
nobs(m1)

## Model 2
m2 <- update(m1, . ~ pct_bypass2_p_t1 + . - bypass_pop_log)

summary(m2)
nobs(m2)


## Collect robust standard errors in list for table
models <- list(m1, m2)
ses <- lapply(models, function(x) sqrt(diag(vcovHC(x, type = "HC1"))))

## Variable names for covariates
covariates <- c("Bypass / Population",
                "Bypass Ratio",
                "Governance Index",
                "Civil Conflicts",
                "Natural Disasters",
                "Years Since Independence",
                "US Ally",
                "No. Major Power Allies",
                "OPEC Member")

## Export to table format 
sink("~/Dropbox/article-manuscripts/baua-isq/baua-isq-final/baua-table-appx-strategic.tex")
stargazer(m1, m2,
          se = ses,
          label = "strategic-models",
          title="Controlling for Strategic Importance of Recipient Country",
          column.labels = c("\\emph{Model 1}", "\\emph{Model 2}"),
          dep.var.caption = c("Unrest Events"),
          dep.var.labels.include = F,
          covariate.labels = covariates,
          omit = c("cowcode", "year", "Constant"),
          p.auto = TRUE,
          t.auto = TRUE,
          digits = 2,
          column.sep.width = c("0pt"),
          notes.append = TRUE,
          notes.align = "l",
          notes.label = "",
          font.size = "footnotesize",
          notes = c("Two-tailed tests. (\\possessivecite{White1980} HC1 standard errors)", 
                    "Negative binomial regression."),
          model.names = FALSE)
sink()






#---------------------------------------------------------------
# Instrumental variables model using bypass per capita measure
#---------------------------------------------------------------

## Model 8, stage 1
mp1 <- lm(bypass_pop_log ~ aid_scandals_t2_major +
            gov_index_t1 + 
            civil_conflict_t1 + 
            nd_count_t1 + 
            mean_global_growth_lag_2 + 
            mean_global_inflation_lag_2 + 
            total_global_disasters_lag_2, 
          data=dat)

## Get data frame for the updated model with extra variables
dat_iv_1 <- model.frame(update(mp1, . ~ . + unrest + year))

## Model 8, stage 2 --- for twostageboot() function
mp2 <- glm.nb(unrest ~ bypass_pop_log +
                gov_index_t1 + 
                civil_conflict_t1 + 
                nd_count_t1  + 
                mean_global_growth_lag_2 + 
                mean_global_inflation_lag_2 +
                total_global_disasters_lag_2,
              data=dat_iv_1, 
              na.action="na.omit")


## Model 9, stage 1
mp1a <- update(mp1, . ~ aid_scandals_t2_extra_continent + . - aid_scandals_t2_major)
summary(mp1a)

dat_iv_2 <- model.frame(update(mp1a, . ~ . + unrest))



## Model 10, stage 1
mp1b <- update(mp1, . ~ scandals_affinity_weighted_past_10 + . + factor(year) - aid_scandals_t2_major - mean_global_growth_lag_2 - mean_global_inflation_lag_2 - total_global_disasters_lag_2)

dat_iv_3 <- model.frame(update(mp1b, . ~ . + unrest + year))

## Model 10, stage 2 --- for twostageboot() function
mp2b <- update(mp2, . ~ . + factor(year) - mean_global_growth_lag_2 - mean_global_inflation_lag_2 - total_global_disasters_lag_2)
summary(mp2b)




#-----------------------------------------------------------
# Uncorrected standard errors models
#-----------------------------------------------------------
## Note: The models above are those that get fed into the twostageboot() function
## to get bootstrapped standard errors.  The second-stage models below are those 
## that appear in the paper. The three chunks immediately following this note
## must be run in succession because they each replace the "predicted_ratio" 
## object for ease of creating a table. The models are the same, but the 
## twostageboot() function automatically performs the transformations below.

predicted_ratio <- predict(update(mp1, data = dat_iv_1))
stage_2_naive_1 <- update(mp2, . ~ predicted_ratio + . - bypass_pop_log)

predicted_ratio <- predict(update(mp1a, data = dat_iv_2))
stage_2_naive_2 <- update(mp2, . ~ predicted_ratio + . - bypass_pop_log)

predicted_ratio <- predict(update(mp1b, data = dat_iv_3))
stage_2_naive_3 <- update(mp2b, . ~ predicted_ratio + . - bypass_pop_log)

## Variable names for table
stage_1_names <- c("Major Aid Scandals",
                   "Extra-Continental Scandals",
                   "Affinity-weighted Scandals",
                   "Governance Index",
                   "Civil Conflicts",
                   "Natural Disasters",
                   "Average Global GDP Growth",
                   "Average Global Inflation",
                   "Total Global Disasters")

## Variable names for table
stage_2_names <- c("Predicted Bypass Per Capita",
                   "Governance Index",
                   "Civil Conflicts",
                   "Natural Disasters",
                   "Average Global GDP Growth",
                   "Average Global Inflation",
                   "Total Global Disasters")

## Get F-statistics on instrumental variables
f1 <- round(summary(mp1)$coefficients[2, "t value"]^2, 2)
f2 <- round(summary(mp1a)$coefficients[2, "t value"]^2, 2)
f3 <- round(summary(mp1b)$coefficients[2, "t value"]^2, 2)

## Export first stage to table format 
sink("~/Dropbox/article-manuscripts/baua-isq/baua-isq-final/baua-table-appx-iv-stage-1.tex")
stargazer(mp1, mp1a, mp1b,
          label = "stage-1-models-appx",
          title = "Stage I - Aid Scandals and Aid Channel Distribution, 2005-2010",
          column.labels = c("\\emph{Model 8, Stage I}", 
                            "\\emph{Model 9, Stage I}", 
                            "\\emph{Model 10, Stage I}"),
          dep.var.caption = c("Logged Bypass Aid Per Capita"),
          dep.var.labels.include = F,
          covariate.labels = stage_1_names,
          omit = c("cowcode", "year", "Constant"),
          p.auto = TRUE,
          t.auto = TRUE,
          digits = 2,
          column.sep.width = c("0pt"),
          add.lines = list(c("Year-fixed Effects", rep("N", 2) ,rep("Y", 1)),
                           c("F-Statistic on Instrument", f1, f2, f3)),
          notes.append = TRUE,
          notes.align = "l",
          notes.label = "",
          font.size = "footnotesize",
          notes = c("Two-tailed tests. Estimated standard errors in parentheses.", 
                    "OLS estimates."),
          model.names = FALSE)
sink()


## Export second stage to table format 
sink("~/Dropbox/article-manuscripts/baua-isq/baua-isq-final/baua-table-appx-iv-stage-2.tex")
stargazer(stage_2_naive_1, stage_2_naive_2, stage_2_naive_3,
          label = "stage-2-models-appx",
          title = "Stage II - Predicted Bypass Per Capita and Unrest, 2005-2010",
          column.labels = c("\\emph{Model 8, Stage II}", 
                            "\\emph{Model 9, Stage II}", 
                            "\\emph{Model 10, Stage II}"),
          dep.var.caption = c("Unrest Events"),
          dep.var.labels.include = F,
          covariate.labels = stage_2_names,
          omit = c("cowcode", "year", "Constant"),
          p.auto = TRUE,
          t.auto = TRUE,
          digits = 2,
          column.sep.width = c("0pt"),
          add.lines = list(c("Year-fixed Effects", rep("N", 2) ,rep("Y", 1))),
          notes.append = TRUE,
          notes.align = "l",
          notes.label = "",
          font.size = "footnotesize",
          notes = c("Two-tailed tests. Estimated standard errors in parentheses.", 
                    "Negative binomial regression models."),
          model.names = FALSE)
sink()



#-----------------------------------------------------------
# Table of bootstrapped coefficients and standard errors
#-----------------------------------------------------------
## Bootstrapped standard errors
set.seed(484)
b1 <- twostageboot(1000, mp1, mp2, "bypass_pop_log", base_data = dat_iv_1, log_t = F)

set.seed(484)
b2 <- twostageboot(1000, mp1a, mp2, "bypass_pop_log", base_data = dat_iv_2, log_t = F)

set.seed(484)
b3 <- twostageboot(1000, mp1b, mp2b, "bypass_pop_log", base_data = dat_iv_3, log_t = F)


## Prepare table
boot_table <- cbind(rep(NA, 3),
                    round(rbind(b1, b2, b3), 3))
colnames(boot_table) <- c(" ","Bypass Coefficient", 
                          "90% lower", "90% upper")
rownames(boot_table) <- c("Major Scandals", "Extra-Continental Scandals", 
                          "Affinity-weighted Scandals")

## Create table
sink("~/Dropbox/article-manuscripts/baua-isq/baua-isq-revised/baua-table-appx-bootstrap-results.tex")
stargazer(boot_table,
          title = "Second-stage Coefficients with Bootstrapped SEs",
          label = "bootstrap-table",
          notes.append = TRUE,
          notes.align = "l",
          notes.label = "",
          notes = c("90\\% CIs from bootstrapped SEs (B = 1000)"),
          digits = 2)
sink()






















